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Despite its apparent simplicity, the idealized model of a particle constrained to move on a circle has 
intriguing dynamic properties and immediate experimental relevance. While a rotor is rather easy to set up 
classically, the quantum regime is harder to realize and investigate. Here we demonstrate that the quantum 
dynamics of quasiparticles in certain classes of nanostructured superconductors can be mapped onto a 
quantum rotor. Furthermore, we provide a straightforward experimental procedure to convert this 
nanoscale superconducting rotor into a regular or inverted quantum pendulum with tunable gravitational 
field, inertia, and drive. We detail how these novel states can be detected via scanning tunneling 
spectroscopy. The proposed experiments will provide insights into quantum dynamics and quantum chaos. 



The classical rotor, a macroscopic particle of mass m confined to a ring, is one of the most studied system in 
classical mechanics. In a gravitational field the rotor becomes a pendulum with its stable equilibrium at the 
bottom. Counter-intuitively, the pendulum can also be stabilized at the top by subjecting the pivot to an 
appropriate drive 1 , in a state referred to as an inverted pendulum. Further interest in classical rotors arose from 
the fact that if one applies periodic kicks to the pendulum, its dynamics becomes chaotic. In contrast to classical 
rotors, quantum rotors are harder to realize 2 . However, with quantum effects coming into play, quantum rotors 
furnish even more abundant phenomena than their classical counterparts. For example, the quantum kicked 
rotor 3 is one of the most featured models of quantum chaos 4 and has various applications 5 " 9 . 

Here we show how a quantum rotor can be obtained in nanostructured superconductors. The idea is to trap 
quasiparticle excitations in a ring-shaped confinement, i.e. a Mexican hat potential, where quasiparticles can form 
Andreev bound states 10 . The desired potential can be achieved in e.g. annular S-N-S junction [see Fig. 1(a)], a 
superconducting film with a ring-shaped surface etch ["blind" ring 11 ; see Fig. 1(b)], or under a perpendicularly 
magnetized ferromagnetic disk [Fig. 1(c)] where screening currents suppress the superconducting order para- 
meter mostly within a ring right under the edge of the disk 1213 . Very recently, it became possible to "draw" the 
suppression of superconductivity practically at will by locally applying electric field in a ferroelectric/super- 
conductor (FE/SC) bilayer via scanning (AFM) tip [see in Fig. 1(d) and Ref. 14 for details]. On the other hand, 
to harbor the Andreev bound states, the ring region should preferably be in the clean limit. As discussed later in 
this manuscript, favorable conditions for experimental detection of the rotor states include: (1) a wide enough 
ring region, in the order of several coherence lengths, and (2) a sufficiently large ring radius when compared to the 
ring width, so that the tunneling through the central potential hill is suppressed. Thus, for conventional super- 
conductors such as Pb or Nb, if a ring is fabricated with its width greater than the e-beam lithography resolution 
(on the order of several tens of nanometers), the ring radius would be on the order of several hundreds of 
nanometers. It would not be trivial to make these conventional superconductors with such long mean free paths, 
although recent experiments by the Roditchev lab 15 show that several monolayers of Pb, when deposited onto a Si 
substrate, become crystalline and exhibit longer mean free paths than bulk Pb. In this respect, graphene- super- 
conductor hybrid is ideal for the realization of our proposal 16 " 18 . Indeed, superconductor can be nanostructured 
with a periodic array of ring perforations (using e.g. e-beam lithography) and covered with a graphene flake [see 
Fig. 1(e)] , so that graphene harbors proximity- induced superconductivity with normal regions above the rings. In 
such a case, the high crystal quality of graphene and therefore the long mean free path (~/mi) would guarantee the 
appearance of the quantum rotor and facilitate its observation. Notably, this would be the very first realization of a 
superconducting quantum rotor in real space, different from the ones realized in Josephson junctions using the 
superconducting U(l) phase space 19 ' 20 . 
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Figure 1 | The 'Mexican hat' order parameter profile required by a 
superconducting rotor state, realizable e.g. in a superconducting (SC) film 
with (a) embedded normal-metal ring (annular S-N-S junction), 
(b) etched "blind" ring on the surface, (c) a ferromagnetic (FM) dot on top 
with out-of-plane magnetization (M), or (d) ferroelectric (FE) film on top 
with ring-shaped domain of inverted polarization (P), or (e) graphene 
flake placed on top of a SC film with ring-shaped perforation. 

Results 

Unusual quasiparticle energy dispersion. In what follows, we focus 
on the case of strongly depleted order parameter on the ring, typical 
for systems of Fig. 1(a), (d) and (e). Nevertheless, our findings can be 
generalized to any other realization of the Mexican hat potential. We 
begin our theoretical analysis of those bound states with a simplified 
model - a step-wise superconducting order parameter with radial 
symmetry A(r), which vanishes between an inner radius R in and 
outer radius R out , but it is otherwise constant and finite. We solve 
the Bogoliubov-de Gennes (BdG) equation in cylindrical coordi- 
nates, and separate the rapidly varying degrees of freedom from 
the slow ones 21 by writing the solutions in the form: 







~fl~ 






-gl- 



H\ l \k P r) + c.c, 



(1) 



where u\ and V\ are BdG wave functions, k F is the Fermi wave vector, r 
is the radial coordinate, is the Hankel function, and//, g\ are 
slowly varying envelope functions. Next, we insert the WKB- 
approximated asymptotic form of the Hankel function and change 
variables 22 23 from r to x, defined by r 2 = x 2 + b 2 . Here b is the WKB 

turning point b= \Jl 2 — 1 /4 j k F , and / is the angular momentum 

quantum number. For large /, the impact parameter is approxi- 
mately b ~ l/kp. We obtain: 



-ivpd x A(x) 
A(x) iv F d x 



(2) 



where A(x) = A(r(x) ) . This maps our two dimensional (2D) problem 
onto a one dimensional (ID) problem. The ID order parameter is 
given by 




if \x\<x in 

if Xj n <\x\ <x out 

if \x\>x out 



(3) 



where x in = yR 2 n — b 2 and x out = y R 2 0Ut — b 2 . In the quasiclassical 

picture, the quasiparticles move along a chord in the x direction. 
The distance b of the chord from the center of the 2D ring 
potential is roughly proportional to the /. This is consistent with 
the classical picture that a particle with higher angular momentum 
is localized further away from the pivot. Each orbit with different 
angular momentum is mapped on a chord of different length. As 
shown in Fig. 2(a), the / ~ 0 state (b ~ 0) resides in the most narrow 
well, with two segments separated by the central hill. As the distance 
b and the angular momentum increases, the chord length of the 
bound state (or, equivalently, the width of the potential well) also 
increases. The bound state energy decreases accordingly as long as 
b < R in . Once b reaches R in , the two segments of the chord merge into 
one wide potential well, causing the energy to drop to a local 
minimum g(/ 0 ), where l 0 ~ k F R in . As b is increased further, the 
length of the chord shrinks and energy increases again. Therefore, 
the dispersion 8(1) must have a local minimum at a non-zero angular 
momentum. This is non-trivial and contradicts classical intuition: for 
such a dispersion faster "rotation" does not always correspond to 
higher energy. Although very different in origin, such a quantum 
effect was found earlier in the roton dispersion in superfluid 4 He 24 " 27 . 

The above reasoning for an unusual dispersion relation is further 
verified by both analytic and numerical solutions using BdG equa- 
tions. For / < l 0 ~ k F R in and A 0 < A^ , the quasiparticle dispersion 
£„(/) can be obtained analytically by matching the solutions at x = 



V F 



2n [ n + - ] + sin 1 (tanh(2^ n )) 



2 {%out ^in ) 



(4) 
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city. For / > / 0 , we obtain 



Here q = -j^-j, £ is the coherence length, and v F is the Fermi velo- 

^X+iY (5) 
2x out \ 2 J 



The energy levels obtained here coincide with Andreev's result 28 for a 
normal slab in contact with two superconductors. Note that s n 
depends on / through x in (l) and x out (l) [see Fig. 2(c)]: it decreases 
with / for / < l 0 and then increases for / > l 0 , as in the case of bound 
states inside a vortex core 22 ' 29 . We confirmed this analytic result 
numerically by solving BdG with B-splines basis (see Methods) 
and direct diagonalization, as shown in Fig. 2(c). 

Quantum rotor. In a realistic system [cf. Fig. 1(c)] one likely finds a 
smooth suppression of the order parameter. Nevertheless, still 
has a local minimum at non-zero angular momentum, as we show in 
Fig. 3(a) [for R = (R out + R in )ll = 12.5^]. The low-lying excitations 
around the dispersion minimum can then be approximated by a 
(I — I ) 2 

quadratic function e(/)~ , where l 0 is the angular 

momentum of the local minimum, and I is the effective inertia. 
The magnitude of I is consistent with m e R 2 , i.e. the inertia 
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Figure 2 | Andreev bound states and their energy dispersion. 

(a) Quasiclassical trajectories of the Andreev bound states in the radially 
step-like order parameter. The bound states are depicted as electrons (solid 
arrow) being Andreev reflected as a hole (dashed arrow). Here b denotes 
the distance of the linear reflection path from the center of the 2D 
potential, (b) Example of a line profile of the order parameter 
corresponding to panel (a). Panel (c) shows the quasiparticle dispersion of 
lowest energy states for the step-like potential featured in panel (b). 

associated to a quasiparticle confined to move along a circular orbit 
with radius R. In the present case, the effective inertia obtained for 
dashed curve in Fig. 3(a) is I = 73 X 10~ 47 kg • m 2 . The effective 
quasiclassical equation of motion, derived from the Hamiltonian 

H= — — — , is very simple, 6 = 0 (6 being the azimuthal angle). 

This equation describes a free rotor - in our case, a quantum rotor. 

To validate the quasiclassical and quadratic approximation, the 
ring is then required to (1) be equipped with a wide enough width 
when compared to the coherence length, so that the energy is low 
enough for those low- lying states, and (2) a large enough radius when 
compared to the width, in order to suppress tunneling across the 
central order parameter hill. Notice that the width of the order para- 
meter well in the ferroelectric/superconductor bilayer can be pre- 
cisely tuned following the method of Ref. 14. This in turn makes 
the shape of e(/) tunable. In other words, one can fine tune externally 
the effective inertia of the quantum rotor we propose here. 

Numerical verification of the rotor state. To confirm these predic- 
tions, we solve numerically the BdG equations on a discrete square 
lattice. Note that in order to compare these results with the 
quasiclassical limit, we have to consider large systems (e.g. 500 X 
500 lattice points). For the size of our nanostructured ring, we take 
the radius R = 100a and width W = 30a, where a = 1 is the lattice 
constant. Here we chose this particular combination of R and W in 
order to have only one peak for / = 0 and a well defined rotor peak for 
/ = Jo- ^ W is larger, multiple overlapping peaks make the analysis 
cumbersome. The value of the superconducting order parameter is 
chosen to be A = O.lf, where t is the hopping amplitude, giving a 
coherence length on the order of the ring width, £ ~ 20a. 

We compare in Fig. 3 the local density of states (LDOS) in the 
middle of planar and ring S-N-S junctions of same width. As seen 
from the graph, the LDOS is modified only within the superconduct- 
ing gap due to multiple Andreev reflections at the normal/super- 



-3-2-10123 

E/A 

Figure 3 | Rotor peak in full BdG calculations. Local density of states 
(LDOS) in the middle of a planar S-N-S junction versus one in the middle 
of a ring junction of same width. The rotor peak at E/A = 0.088 is present 
only in the ring junction. 

conducting interface. For both situations there is a main peak 
around E/A = 0.55 which is mainly given by paths of length W 
and corresponds to the state with / = 0. Interestingly, for the ring 
geometry an additional peak appears at low energies. This peak is the 
low-lying state with / = l 0 discussed in previous section and is given 
by the longest chords in the ring. In other words, this peak corre- 
sponds exactly to the rotor state. 

Discussion 

In this Section, we present further implications of the existence of the 
rotor state, and its manifestations in the presence of applied drive. 
These effects are interesting in their own right, but also essential for 
the experimental verification of the rotor state. 

Quantum pendulum. It is well-known that superconducting pro- 
perties can be strongly affected by an electric current. An applied dc 
current in our system leads to the following effective Hamiltonian for 
the low-lying states: 

H=^+v,k, (6) 

where k = (k F cos 0, k F sin 9) and v s is the superfluid velocity 
corresponding to the applied transport current. The equation of 
motion then exhibits nonlinear dynamics: 

kpVs . 
I 



0- 



-sin0 = O. 



For our choice of coordinates, this corresponds to an inverted 
pendulum. Nevertheless, the "gravitational" force drives the 
system away from the 9 = 0 (up) unstable equilibrium point. To 
analyze the behavior of this quantum pendulum, we use the effective 
Hamiltonian in Eq. (6), and show that the thermally averaged zero 

E„l«A„(r)| 2 <5(£-£„) rf£) 



energy local density of states, 5V « ; — . 

J _oo 4Tcosh (e/2T) 

has a peak at the bottom position, 0 = n [see Fig. 4(c)]. Hence in 

the presence of an applied dc current our system becomes the 

quantum analogue of a pendulum in a gravitational field. 

Numerical results with applied current. We next show the full 
numerical solutions of the BdG equations on a discrete lattice in 
the presence of an external dc current. A supercurrent is applied to 
the system in y- direction by imposing a phase gradient, d^/dy, on the 
superconducting order parameter. 
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Figure 4 | Manifestation of the quantum rotor in the local density of states, (a) Parabolic approximation (dotted curve) of the energy dispersion [for 
realistically smoothened order parameter profile shown as inset] . The resulting Hamilton equation of motion corresponds to a rotor, (b) The local density 
of states (LDOS) calculated for a quantum rotor, for the same parameters as in (a), (c) The calculated LDOS for a quantum pendulum, obtained by 
injecting a dc current in the indicated direction ( v s = 1.0 m/s). The finite transport current is effectively equivalent to gravity in the classical pendulum 
picture. 



The dashed line in Fig. 5(a) shows the calculated LDOS(E) for the 
case without applied current, taken at the bottom position of the ring 
(at 0 = n). One notices two main peaks, largest of which corresponds 
to the state with / = 0, whereas the smaller one is the state with 1= l 0 . 
The solid curve in Fig. 5(a) is calculated at the same place on the ring, 
but now in the presence of the current along the ^-direction. 
According to Fig. 4, a pendulum state is expected. We indeed observe 
that the rotor peak is split into pendulum states since there is a phase 
difference along the chords which contribute to those states. 
However, because the chiral symmetry is not broken, states with 
opposite chirality are located at top and bottom locations respect- 
ively, as shown in the LDOS map in Fig. 5(b) [for energy correspond- 
ing to the lowest energy peak in Fig. 5(a)]. In other words, the two 
pendulum states at 0 = n and <j> = 0 are degenerate in energy, but the 
application of a small perpendicular magnetic field breaks the chiral 
symmetry and selects only one of them [as shown in Fig. 5(c)]. The 
LDOS maps for all the energies and other field/current configura- 
tions can be found in the Supplementary Material. 

Inverted pendulum. Finally we return to Eq. (6) to point out that 
although the inverted pendulum is unstable, it can be driven to 
stability by horizontal or vertical pivot motion (cart-pole or 
Kapitza pendulum, respectively 30 ). If the pivot is vibrated with an 
amplitude A cos cot, and the frequency co is much higher than the 

pendulum's natural frequency co N , i.e. more precisely co^> — co N , 

where w N = \/hk F v s /I, the effect of the fast pivot motion will be 

identical to an effective potential 30,31 V e jf = - IR 2 a 2 co 2 sin 2 0. For 

such a case, which in a system of Fig. 1(c) [(d)] can be realized by 
an ac in-plane magnetic [electric] field, we calculated the local 
density of states [see Fig. 6] and indeed observed a peak at the top 
position 6 = 0 of the pendulum, corresponding to the inverted 
pendulum state. Note, however, that the peak at 6 = n for a 
regular pendulum state remains, which is different from the 
classical case. Quantum mechanics allows the particle to reside in 
both potential wells, while classical mechanics forces the particle to 
choose either bottom or top localization. Finally we note that while 
the inverted pendulum state is stabilized by a high frequency drive, 
low- frequency drives are also of interest as they lead to the classically 
chaotic regime. Therefore, the quantum rotor proposed here is a 
testbed for further studies of quantum chaology 32 . Applying pulsed 
current to our device would transform the superconducting rotor 
into a pendulum with pulsed gravity, i.e. another quantum chaotic 
system 7 " 9 . 

In summary, we provided analytical and numerical results for 
novel quantum rotor states in nanostructured superconductors. 
Besides being a remarkable example of a quantum analogue of a 



classical system, the superconducting rotor (i) can be realized in 
various superconducting systems, (ii) has a tunable inertia and grav- 
itational field, (iii) can be externally manipulated through effective 
tilt, pulsed gravity, and pivot oscillations, and (iv) can be converted to 
a quantum pendulum, inverted pendulum, or be driven to a chaotic 
regime. Hence this quantum rotor has the potential to provide 
insights into a variety of phenomena which certainly deserve further 
experimental and theoretical investigation. 

Methods 

The Bogoliubov-de Gennes (BdG) equations are used to calculate the spectrum and 
LDOS throughout this work. The BdG equations provide mean field description of 
superconductivity and are generalized from the Hartree-Fock equations. We first 
applied the quasiclassical approximation, which separates the fast oscillation mode 
with length scale l/k F . It allowed us to discuss physical quantities varying over the 
superconducting coherence length, and also provided clear physical picture of the 
quasiclassical trajectories. For the quasiclassical limit, the numerical calculations are 
done within the B-spline basis 33 " 35 , which is very efficient for inhomogeneous pro- 
blems such as the one presented here. 

We also performed full BdG calculations for a discrete square lattice by using the 
Chebyshev-BdG method 18,36 . This method can extract relevant physical quantities, for 
example LDOS, from large systems (here we used a square lattice with 500 X 500 grid 
points) and avoid large matrix diagonalization. 

In order to describe inhomogeneous s-wave superconductivity we use a real-space 
tight-binding formulation. The corresponding Hamiltonian written in Nambu spinor 
form is the following: 



2.5 
2 



8 15 
3 1 

0.5 
0 



(a) ' 


v s *b — - 














Figure 5 | Realization of the quantum pendulum, (a) Full BdG calculation 
of the LDOS vs. energy in 6 = % location on the ring (of width W = f and 
radius R = 5£), for the quantum rotor (no current) and pendulum state 
(with current generated via d(f>/dy = 0.0033/a). The 2D LDOS for a 
quantum pendulum is shown for El A = —0.088 in the absence (b) and 
presence (c) of a small perpendicular magnetic field, O/O 0 = 0.15. 
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Figure 6 | Pivot oscillations and the inverted pendulum. The 

localization of a regular and an inverted pendulum, directly 
observable in STM imaging of LDOS, in the case of the oscillating pivot 

with amplitude of the oscillation A = 10 nm, for high drive frequency 
cd = 50 GHz^> — co w ~l GHz. The inset shows LDOS along the circular 
path of the rotor. The blue/red curve corresponds to the LDOS without/ 
with ac drive. 



where is a 2 X 2 matrix: 



^'4 > + w >-*>■ 



where the sum is over nearest and next-nearest neighbours. The superconducting 
order parameter, A,-, is of s-wave spin-singlet type and is finite outside the ring and 
and zero inside it. The hopping parameters for nearest neighbor, t = 2/3 1 0 , and next 
nearest neighbor, t' = V6t 0 together with the chemical potential /u= —2t Q are chosen 
such that the Fermi surface is circular and therefore there is no anisotropy in the 
Fermi velocity. t 0 is an arbitrary energy unit which is set to 1. 

We considered both the influence of an in-plane current and/or an out-of-plane 
magnetic field on the rotor states. The supercurrent is described by a phase gradient 
imposed on the superconducting order parameter. The out-of-plane magnetic field is 
considered to be weak, the flux through a circle of radius R was typically taken to be O 
= 0. 15® 0 , so that Meissner currents are negligible and a self- consistent solution is not 
necessary. The magnetic field modifies the hopping parameters in the usual way 

r e v -^-> ' 

through the Peierls phase factors, ty = \ ty\exp i— Adl 

For finding the local density of states inside the ring we use the Chebyshev-BdG 
method which provides an efficient way of calculating the Gorkov Green function: 



Gy(ca) = {vac\^ ^ jG(a))(c^ c^^vac), 



(10) 



where ±G(co + irj) = [co + irj — H] 1 , | vac) is the vacuum and the diagonal( normal) 
and off-diagonal(anomalous) components can be expressed as: 



Gj>) = (c it |G(a>)c+ 



5J» = 



4 Igko)) 4, 



(ii) 



(12) 



Each of these components can be expanded in terms of Chebyshev polynomials 
and the moments of the expansion can be efficiently calculated through a recursive 
procedure: 



=,11(12) 



(&)- 



where the coefficients a l n l (i,j) = (c^ | v n ) and ±a l „(i,j) = (c]± | v„^ can be obtained 

by an iterative procedure involving repeated applications of the Hamiltonian on 
iterative vectors \v n ): 

\v n+ i)=2H\v n )-\v n - l ), (14) 

where |vo) = |cj^ and |v_i) = 0. The local density of states (LDOS) is then: 

1 



Ni{co) = --Im[G]l(co)]. 



(15) 



We performed the expansions separately for each grid point and further accelerated 
our calculations by using graphical processing units (GPUs) to perform the sparse 
matrix - vector multiplications. 
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